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(54) Abstract Title 

Method of extracting a signal from a contaminated signal 

(57) A method is provided of extracting desired signals s^ from contaminated signals measured via 
respective communication channels. The system comprising the desired signals and the channels is 
modelled as a state space model, fn the model, the desired signals have time-varying characteristics which 
vary more quickly than second time-varying characteristics of the channels. The method may be employed to 
extract individual speech signals from speakers s^*^* and St*^* in a room 1, contaminated signals y^ being 
sampled by microphones 6, 7, 8 and processed by a computer 9 which extracts the individual speech signals 
for supply to respective loud speakers 13, 14. 
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At least one drawing originally filed was informal and the print reproduced here is taken from a later filed formal copy. 
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METHOD OF EXTRACTING A SIGNAL 

The present invention relates to a method of extracting a signal. Such a method may be 
used to extract one or more desired signals from one or more contaminated signals 
received via respective communications channels. Signals may be contaminated with 
noise, with delayed versions of themselves in the case of multi-path propagation, with 
other signals which may or may not also be desired signals, or with combinations of 
these. 

The communication path or paths may take any form, such as via cables, 
electromagnetic propagation and acoustic propagation. Also, the desired signals may in 
principle be of any form. One particular application of this method is to a system in 
which it is desired to extract a sound signal such as speech from contaminating signals 
such as noise or other sound signals, which are propagated acoustically. Another 
particular application of this method is to digital communication, for example where the 
signals represent digitised data and the or at least one channel is time-varying. 

According to a first aspect of the invention, there is provided a method of extracting at 
least one desired signal from a system comprising at least one measured contaminated 
signal and at least one communication channel via which the at least one contaminated 
signal is measured, comprising modelling the system as a state space model in which 
the at least one desired signal has first characteristics and the at least one 
communication channel has second characteristics which are different from the first 
characteristics. 

State space models are known in mathematics and have been applied to the solution of 
some practical problems. A state space model relates to a system in which there is an 
underlying state of the system which it is desired to estimate or extract. The state is 
assumed to be generated as a known function of the previous state value and a random 
error or disturbance term. The available measurements are also assumed to be a known 
function of the current state and another random error or noise term. 
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It has been surprisingly found that a state space model may be successfully applied to 
the problem of extracting one or more desired signals from a system comprising one or 
more measured contaminated signals and communication channels. It has been realised 
that the time varying characteristics of the or each desired signal differ from the time- 
varying characteristics of the or each communication channel in such a way that the or 
each desired signal can be extracted from the or each contaminated signal which is 
actually measured. This technique makes possible the extraction of one or more desired 
signals in a tractable way and in real time or on-line. A further advantage of this 
technique is that future samples are not needed in order to extract the samples of the 
desired signal or signals although, in some embodiments, there may be an advantage in 
using a limited selection of future samples. In this latter case, the samples of the desired 
signal or signals are delayed somewhat but are still available at at least the sampling rate 
of the measured signals and without requiring very large amounts of memory. 

The first characteristics may be fixed and known. For example, this techniques may be 
applied to echo cancellation. As an alternative, the first characteristics may be lime- 
varying. Although the second characteristics may be unknown, they may be fixed and 
known. For example, the characteristics may have been previously determined by 
training data or in any other appropriate way. 

The second characteristics may be time-varying. The first time-varying characteristics 
may vary on average more quickly than the second time-varying characteristics. For 
many systems, the characteristics of the desired signal or signals vary relatively rapidly 
whereas the characteristics of the communication charmel or channels vary more slowly. 
Although there may be abrupt changes in the channel characteristics, such changes are 
relatively infrequent whereas signals such as speech have characteristics which vary 
relatively rapidly. By modelling these characteristics in such a way that the different 
rates of variation are modelled, the extraction of one or more signals is facihtated. 

The at least one communication charmel may comprise a plurality of conununication 
charmels and the at least one contaminated signal may comprise a plurality of 
contaminated signals. The at least one desired signal may comprise a plurality of 
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desired signals. The number of communication channels may be greater than or equal 
to the number of desired signals. Although it is not necessary, it is generally preferred 
for the number of measured signals to be greater than or equal to the number of desired 
signals to be extracted. This improves the effectiveness with which the method can 
recreate the desired signal and, in particular, the accuracy of reconstruction or extraction 
of the desired signal or signals. 

The at least one contaminated signal may comprise a linear combination of time- 
delayed versions of at least some of the desired signals. The method is thus enable of 
extracting a desired signal in the case of multi-path propagation, signals contaminating 
each other, and combinations of these effects. 

The at least one contaminated signal may comprise the at least one desired signal 
contaminated with noise. Thus, the method can extract the or each desired signal from 
noise. The at least one channel may comprise a plurality of signal propagation paths of 
different lengths. 

The at least one desired signal may comprise a sound signal. The at least one sound 
signal may comprise speech. The contaminated signals may be measured by spatially 
sampling a sound field. For example, acousto-electric transducers such as microphones 
may be spatially distributed in, for example, a room or other space and the output 
signals may be processed by the method in order to extract or separate speech from one 
source in the presence of background noise or signals, such as other sources of speech 
or sources of other information-bearing sound. 

The at least one desired signal may be modelled as a time-varying autoregression. This 
type of modelling is suitable for many types of desired signal and is particularly suitable 
for extracting speech. As an alternative, the at least one desired signal may be 
modelled as a moving average model. As a fiirther alternative, the at least one desired 
signal may be modelled as a non-linear time- varying model. 
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The at least one communication channel may be modelled as a time-varying finite 
impulse response model. This type of model is suitable for modelling a variety of 
propagation systems. As an alternative, the at least one communication channel may 
be modelled as an infinite impulse response model. As a further alternative, the at least 
one communication channel may be modelled as a non-linear time- varying model. 

The state space model may have at least one parameter which is modelled using a 
probability model. The at least one desired signal may be extracted by a Bayesian 
inversion. The Bayesian inversion may be performed by a sequential Monte Carlo 
method. These techniques are particularly effective for complex models and are 
potentially implementable on parallel computers. 

According to a second aspect of the invention, there is provided a program for 
controlling a computer to perform a method as claimed in any one of the preceding 
claims. 

According to a third aspect of the invention, there is provided a carrier containing a 
program in accordance with the second aspect of the invention. 

According to a fourth aspect of the invention, there is provided a computer programmed 
by a program according to the second aspect of the invention. 

The invention will be further described, by way of example, with reference to the 
accompanying drawings, in which: 

Figure 1 is a diagram illustrating a signal source and a communication channel; 

Figure 2 is a diagram illustrating audio sound sources in a room and an apparatus for 
performing a method constituting an embodiment of the invention; and 

Figure 3 is a flow diagram illustrating a method constituting an embodiment of the 
invention. 
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In order to extract a desired signal, a parametric approach is used in which the data are 
assumed to be generated by an underlying unobserved process described at time t with 
the variable Xt. The variable x, contains information concerning the waveforms of the 
different sources. The problem of extracting a desired signal may be too complex for a 
good deterministic model to be available, as little is known concerning the real structure 
of the problem. Altematively, if a good deterministic model is available, this may lead 
to a large set of intractable equations. 

The models for the way sound is generated give expressions for the likely distribution 
of the current *state' Xt given the value of the state at the previous time step Xt-i. This 
probability distribution (known as the state transition distribution, or *prior') is written 
as p(x,|Xt.i). How the current observed data depend on the current state is specified 
through another probability distribution p(yi|xt) (known as the observation distribution, 
or likelihood'). Finally, how the state is likely to be distributed at the initial time 
instant is specified by p(xo). Specific forms for these three distributions are given later. 

A solution to circumvent the problems mentioned earlier comprises introducing 
uncertainty in the equations through probability distributions. More precisely this 
means that, instead of assuming, in a discrete time set-up, that Xt+i is a deterministic 
function of past values e.g. Xt+i = Axt where A is a linear operator, the plausible regions 
of the state space where the parameter can lie are described with a conditional 
probability distribution p (xt+i = % \ x,), the probability of x^+i being equal to x given the 
previous value Xt. This may be expressed, for example, as Xt+i = Axt + Vt, where Vi is an 
uncertainty or error distributed according to a Gaussian distribution around zero. The 
way the distribution is spread indicates the degree of confidence in the deterministic 
component Axt; the narrower the spread, the greater the confidence. This type of 
modelling proves to be robust in practice to describe very complex processes while 
being simple enough to be used in practice. 

As mentioned above, whereas the structure of the process is assumed known (because 
of plausible physical assumptions), the variable Xi is not observed and solely 
observations yt are available. In general, the observation mechanism, i.e. the 



transformations that the process of interest undergoes before being observed, will 
depend on the variable Xt, but again some randomness needs to be incorporated in the 
description of the phenomenon, for example to take into account observation noise. 
Again this is done by means of a probability distribution p (yt=y|xt), the probability of yt 
being equal to y given that the parameters of the underlying process are represented by 
Xt- In the previous example, a possible observation process is of the form yt=Cxt+Wt, 
that is some transformation of the parameter Xt describing the internal state of the 
process corrupted by an additive observation noise Wt- 

In the case of audio source separation, the parameter Xt contains the value St of the 
desired waveforms of the sources at time t and the evolving parameters of the sources 
(a,) and the mixing system (ht). This is illustrated diagrammatically in Figure 1 for a 
single communication chaimel. Transition probabihty distributions are defined by: 

p(at+i = alat) 
p(ht+i = h|ht) 

and as it is not known how to characterize this evolution, apart from the fact that at is 
expected to evolve rapidly compared to ht (which might also evolve abruptly, but this is 
expected to happen rarely): 

where the distributions of , vf are each isotropic, but have a different spread to 
reflect the different non-stationarity time scales. The set of the parameters at may be 
thought of as evolving spectral characteristics of the sources and the parameters ht as 
classical finite impulse response (FIR) filters. Now that the parameters of the system 
are described, the waveforms may be assumed to evolve according to a probability 
distribution: 
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p(st+i = s|st, at) 
which describes the model of the sources. 

One example of audio signal separation or extraction is where the sources are speech 
sources. The parameter Xt may then be modelled as a time-varying autoregressive 
process, for example by the expression: 

p ■ - . 

i=t 

where the parameters aj^t are filter coefficients of an all-resonant filter system 
representing the human vocal tract and Ct represents noise produced by the vocal cords. 

Figure 2 illustrates a typical application of the present method in which two desired 
sound sources, such as individuals who are speaking, are represented by s^'\ and s^^\ 
located in an enclosed space in the form of a room 1 having boundaries in the form of 
walls 2 to 5. Microphones 6 to 8 are located at fixed respective positions in the room 1 
and sample the sound field within the room to produce mesisured contaminated signals 

y^^\y y^^V The measured signals are supplied to a computer 9 controlled by a 
program to perform the extraction method. The program is carried by a program carrier 
10 which, in the example illustrated, is in the form of a memory for controlling the 
operation of the computer 9. The computer extracts the individual speech signals and is 
illustrated as supplying these via separate channels comprising amplifiers 1 1 and 12 and 
loudspeakers 13 and 14. Alternatively or additionally, the extracted signals may be 
stored or subjected to further processing. 

Figure 2 illustrates some of the propagation paths fi-om the sound sources to one of the 
microphones 7. In each case, there is a direct path 15, 16 ft^om each sound source to the 
microphone 7. In addition, there are reflected paths firom each sound source to the 
microphone 7. For example, one reflected path comprises a direct sound ray 1 7 fr^om 
the source which is reflected by the wall 2 to form a reflected ray 18. In general, there 
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are many reflected paths from each sound source to each microphone with the paths 
being of different propagation lengths and of different soimd attenuations. 

The objective is to estimate the sources Si given the observations yt and, to be of 
practical use, it is preferable for the estimation to be performed on line with as little 
delay as possible. Whereas the modelling process describes the causes (x,) that will 
result in effects (the observations yt), the estimation task comprises inverting the causes 
and the effects. In the frameworic of probability modelling, this inversion can be 
consistently performed using Bayes' theorem. In the context of statistical filtering, that 
is the computation of the probability of Xt given all the observations up to time t, namely 
PWyi:t), given the evolution probability p(xi|xt.,) and the observation probability 
P(yt|xt), the application of Bayes' rule yields 



p(^r \yuo ) = -—I — 



This equation is fundamental as it gives the recursion between the filtering density at 
time t-1, i.e. p(Xt-i|yi:,-i), and the filtering density at time t, p(xt|yi:t). The problem is 
that, in practice, the integral cannot be computed in closed-form and it is desired to 
compute these quantities in an on-line or real time manner. 



Whereas the use of probability distributions can be viewed as a way of representing in a 
parsimonious way the concentration of members of a population in certain regions of 
their feature space, Monte Carlo methods work in the opposite direction and rely on the 
idea that a probability distribution can be represented with an artificially generated set 
of samples distributed according to this distribution. This requires that the 
concentration of the samples in a given zone of the space of features is assumed to be 
representative of the probability of this zone under the distribution of interest. As 
expected, the larger the population, the more accurate the representation is. This 
approach possesses the advantage in many cases of greatly simplifying computation, 
which can to a large extent be performed in parallel. 
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The following represents a basic explanation of the techniques involved in the present 
method. This is followed by a detailed desaiption of a specific embodiment. 

It is assumed that, at time t-1, there are N » 1 members of a population of "particles" 
x,^l{ distributed according to the distribution p(xi.i|yi:i-i). Then it is possible to guess 
where the particles are going to evolve using the evolution distribution p(xt|xt.i) or prior, 
from which it is typically easy to sample. Thus, * scouts' are being sent into the regions 
which are likely, according to the prior on the evolution of the parameters. When the 
next observation yi is available, the prediction needs to be corrected as the 'scouts* did 
not take into account the information brought by the new observation which can be 
quantifyied with p(yt|xi). Some of the particles will be in regions of interest but 
typically insufficient numbers, or there might also be too many members of the 
population in regions where many less should be. 

A way of regulating the population consists of multiplying members in underpopulated 
regions and suppressing members in overpopulated regions. This process is refered to 
as a selection step. It can be proved mathematically that the quantity that will decide 
the future of a 'scout' is the importance function: 

and valid mechanisms include taking the nearest integer number to Nw^*^ to determine 

the number of "children" of scout number i or randomly choosing the particle x}'^ with 

probability w^'^ . After the selection process, the children are approximately distributed 

according to p (xilyin) and the next data can be processed. This is a very general 
description of the algorithm, and many improvements are possible as described 
hereinafter for a very specific case. 
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A general recipe for how to do sequential Monte Carlo estimation is given. This is the 
most basic form of the method and is described in NJ. Gordon, DJ. Salmond and 
A.F.M. Smith, "Novel approach to nonlinear/non-Gaussian Bayesian state estimation", 
lEE Proceedings-F, vol. 140, no. 2, pp. 107-1 13, 1993, the contents of which are 
incorporated herein by reference. 

In an initial step 20, time t is set to zero. N so-called "particles" are randomly chosen 
from the initial distribution p(xo) and labelled x^^'^o, where k is an integer from 1 to N. 
The number N of particles is typically very large, for example of the order of 1,000. 
The steps 20 and 21 represent initialisation of the procedure. 

State propagation is performed in a step 22. In this step, for each particle from the 
current time t=l , an update to the cuirent time is randomly chosen according to the 
distribution p(xt+i|x^*'\). The updated particles x^\+\ have plausible values which lie in 
the expected region of the space according to the state transition distribution. 

In a step 23, the next measured value y^+i is obtained from measuring or sampling the 
sound field. In a step 24, the weights of the particles are calculated. Because the 
updated particles x^^Vi were generated without reference to the new measurement yi+i. it 
is necessary to calculate a weight between zero and one for each particle representing 
how "good" each particle actually is in the light of the new measurement The correct 
weight value for each particle is proportional to the value of the observation distribution 
for the particle. However, the sum of all of the weights is required to be equal to one. 
Accordingly, the weight of each particle x^^\+] is given by: 

A step 25 then reselects or resamples the particles in order to return to a set of particles 
without attached weights. The step 25 reselects N particles according to their weights 
as described in more detail hereinafter. A step 26 then increments time and control 
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returns to the step 22. Thus, the steps 22 to 26 are repeated for as long as the method is 
in use. 

The value of each sample may, for example, be calculated in accordance with the 
posterior mean estimator as: 

The techniques described hereinbefore are relatively general and are not particularly 
limited to the case of separating or extracting speech or other audio signals from 
"contaminated" measurements. A method directed more explicitly to extracting audio 
source signals will now be described in greater detail. 

The sources are assumed to be modelled with time varying autoregressive processes, i.e. 
source number i at time t is a linear combination of pi (the so-called order of the model) 
past values of the same source (Si,t-i,...^Si,t-pi, or in short and vector form s,>i:t-pi) 
perturbed by a noise, here assumed Gaussian, vf , : 

This type of process is very flexible and allows for resonances to be modelled. The 
coefficients aj^tof the linear combination, are time- varying in order to take into account 
the non-stationarity of the resonances of speech: 

(2) 



The observations, which consist of the superimposition of the different sources 
(possibly delayed) at microphone j at time t are assumed generated by the following 
process 
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(3) 



i.e. the observation is the sum over all the sources of filtered versions of the sources (the 
filter is of length lij fi"om source i to microphone j, and introduces delays) perturbed by 
an observation noise wj^t. The transfer fiinction firom source i to microphone] evolves 
in time according to the equation 



It is assumed that the characteristics of v*,^, and wj^t are known, whereas in the full 
version of the algorithm these parameters are also estimated. 

The steps in the sequential Monte Carlo algorithm for audio source separation are 
essentially as shown in Figure 3 but are changed or augmented as described hereinafter. 
The state vector xi is defined as a stacked vector containing all the sources, 
autoregressive parameters and transfer functions between sources and microphones. 

The initial values are randomly chosen from a normal distribution with a large spread 

In the step 22, random noise v^*^\t+i, v*\,+i and v^^^ij^i+i are generated from Gaussian 
distributions, for example as disclosed in B.D. Ripley, "Stochastic Simulation'*, Wiley, 
N.Y. 1987, the contents of which are incorporated herein by reference. These quantities 
are then added to their respective state variables: 




CI, 



r./ + I 
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The step 24 calculates un-normalised weights as follows: 



exp 



• • X exp 



-1 

2^1 



'n.t+l 



<*) 
i.lJ-lj.„*\ 



The weights are then renormalised in accordance with the expression: 



The variability of the variances of the excitation noise and the observation noise may be 
taken into account by defining evolution equations on the <|>;^+,.^ = log (a and 

<l>y.,*..v = log («J as follows 



where vj';^, and v^%, are i.i,d, sequences distributed according to Gaussian 
distributions. 

It is then possible to improve the performance with some or all of several modifications 
to the basic algorithm as follows: 

(a) Estimation of the time- varying noise variances. Generate random noise 
>^iv+i from normal distributions and compute 
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which are transformed with an exponential to obtain the variances of the noises, i.e 



(b) By using the structure of the model, the mixing filter coefficients and the 
autoregressive coefficients can be integrated out. The implication is that the calculation 
of the weight is modified to 



This weight is calculated sequentially by use of a standard procedure, the Kalman filter, 
(see Appendix Eq. (50)) applied to the state space model for which the state consists of 
the autoregressive coefficients and mixing filter coefficients shown in the following 
equations (1 7) to (19). The advantage of this method is that the number of parameters 
to be estimated is significantly reduced and the statistical efficiency is much improved. 

(c) The distribution that propagates the values of the sources is modified to an 
approximation of the filtering density p{c^\yt:^a^,^\^^ which, in 

contrast to the basic algorithm, takes into account the new observation yt, hence leading 
to improved statistical efficiency once again. This density is the byproduct of a Kalman 
filter (say Kalman filter #2) applied to a state space model for which the state consists of 
the sources Xt and the parameters shown in the following equation (14) depend upon the 
filtered mean estimate of the mixing coefficients h and autoregressive coefficients a, 
previously obtained from Kalman filter # 1 (see Appendix) 

(d) Diversity among the particles is introduced by using a Metropolis-Hastings update 
on the particle whose target distribution is pid{x.al^,^^ , cr^ , ^ \y^.^ ) 

This step is introduced after the step 25 and the details are given hereinafter. 



= exp(^„.^,^,)and ct'^ ,,, =exp 
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(e) Estimation is delayed for improvement purposes, i.e. to compute the value of the 
sources at time t, wait and take into account some future observations yr+y...., yt+t for 
some L, Such a technique is called fixed-lag smoothing and does not modify the 
algorithm as it is merely necessary to wait for the data yt+i 

The full details of these different steps are given hereinafter. 



The whole system including sources and channels is modelled as a state space system, a 
standard concept from control and signal processing theory. The definition of a state 
space system is that there is some underlying state in the system at time ^ denoted Xi, 
which it is desired to estimate or extract. In the present case, this comprises the 
underlying desired sources themselves. The state is assumed to be generated as a 
known function of the previous state value and a random error or disturbance term: 



where At+i(. . .) is a known function which represents the assumed dynamics of the state 
over time. Similarly the measurements at time t, denoted yt, are assumed to be a known 
function of the current state and another random error or noise term: 



where (...) is a known function representing the contaminating process. In the 
present case, (. . .) represents the filtering effects of the channel(s) on the source(s) 
and represents the measurement noise in the system. 

A key element of the present method is that both the source(s) and channel(s) have 
different time-varying characteristics which also have to be estimated. In other words, 
the functions and themselves depend upon some imknown parameters, say 
9/ and 0 f . In the present case, Q^* represents the unknown time- varying autoregressive 
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parameters of the sources and G , represents the unknown time-varying finite impulse 

response filter(s) of the chaimel(s). These time varying characteristics are modelled 
with additional state transition functions, as described below. 

The problem addressed is the problem of source separation, the n sources being 
modelled as autoregressive (AR) processes, from which, at each time t, there are m 
observations which are convolutive mixtures of the n sources. 
Source / can be modelled for /=1, . . . as: 



and Pi is the order of the AR model. It is assumed that si^-pi +i:o - 

N{ml , ) for / = 1, n. (vf, is a zero mean normalized i.i.d. Gaussian 

sequence, i.e. for / = 1 . and / = 1 , . . . 



The (cj ^ . , )i ^ 1,. .,n are the variances of the dynamic noise for each source at time t. 

It is assumed that the evolving autoregressive model follows a linear Gaussian state- 
space representation. 




(4) 




-N(0,1) 



(5) 



4»v.,V*l 



(6) 





(7) 
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and ai.o -N{m% ) . The model may include switching between different sets of 
matrices ^4/,^/ to take into account silences, for example. Typically A? =/^. and 



The mixing model is assumed to be a multidimensional time varying FIR filter. It is 
assumed that the sources are mixed in the following manner, and corrupted by an 
additive Gaussian i.i.d. noise sequence. At the yth sensor, and fory = 1 , „ ., m: 

n 

= S Kjjhju-iu^l + <^ w,,V ^yV (8) 

where Ifj is the length of the filter from source / to sensor j. The series (uy^,=.y ... is a 
zero mean normalized i.i.d. Gaussian sequence, i.e. fory = 1 . n and t = 1 , . 



M0,1) (9) 



The quantities (a Ij, ),^^ ^ are the variances of the observation noise for each sensor at 

time L The Wy.,are assumed independent of the excitations of the AR models. The 
constraint [h/j.,]i,i = 1 and [h/y.Jk, i = 0 fory = / mod /w and z = 1,..., n are imposed. In 
the case m = n, this constraint corresponds to [hij,(]i.i = 1 and [hij.t]k,i = 0 for j = L As 
for the model of the sources, it is assumed that the observation system also follows a 
state-space representation. Writing ^^jj ^log (a ^ ^> , ) : 



- J* 



with 



<j-' "~^Kj..AjUr, Mo,i) (11) 



and \,,o~n{<Po] 
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For each of the sources /, the signal can be rewritten in the following form: 



where A, A max max {?,. ^.} and A*, is a x Xj matrix defined 



as: 



(12) 



The dynamic system equations can be rewritten as: 



with 



V,' = (vm . - , y , x,2 (5,';, „ . . . 

[c'j contains the mixing system 

and w' =(w,,...w„ J ,D';-diag(a^,...a^„). 



(13) 



(14) 



(15) 



Defining the "stacked" parameter vectors 



[A,]y-yA-'/ +;i4A[\^.,[,/ = l,...,"»;7, =l,..,";^.= l....,/(, (16) 



and AX!"*! /^i* S>-i 2,-1 'v ~ S>-i i possible to consider the following state 
space representations. 
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y,=C^h, +D^w^ 
where A^B^ e , C* g i?""'* , D,* e 



(17) 



(18) 



and 



X, =C>, +/)>; (19) 



with € /?''"» C; ^R'"'-where 



(20) 



This is a practical interest as it allows the mixing filters and autoregressive filters to be 
integrated out because, conditional upon the variances and states x„ the systems (17) 
and (19) are linear Gaussian state space models. Together with (14) they define a 
bilinear Gaussian process. 

Given the number of sources m, pi and ly it is required to estimate sequentially the 
sources (x,),=,,. . . and their parameters G aIo A ct ^ \ 

from the observations More precisely, in the framework of Bayesian estimation, 
one is interested in the recursive, in time, estimation of posterior distributions of the 
type 

p(dQ , dx,\y^,^^ ): when Z = 0 this corresponds to the filtering distribution and when 
Z, > 0 this corresponds to the fixed-lag smoothing distribution. This is a very complex 
problem that does not admit any analytical solution, and it is necessary to resort to 
numerical methods. Such a numerical method is based on Monte Carlo simulation. 
Subsequently the following notation will be used: 
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A simulation-based optimal filter/fixed-lag smoother is used to obtain filtered/fixed-lag 
smoothed estimates of the unobserved sources and their parameters of the type 

/.(/;)Aj/,(e„x,)/7(rfe„.&,|y„,J (21) 

The standard Bayesian importance sampling method is first described, and then it is 
shown how it is possible to take advantage of the analytical structure of the model by 
integrating onto the parameters a^ and ht which can be high dimensional, using Kalman 
filter related algorithms. This leads to an elegant and efficient algorithm for which the 
only tracked parameters are the sources and the noise variances. Then a sequential 
version of Bayesian importance sampling for optimal filtering is presented, and it is 
shown why it is necessary to introduce selection as well as diversity in the process. 
Finally, a Monte Carlo filter/fixed-lag smoother is described. 

For any ft it will subsequently be assumed that |/^ (/, \< + oo. Suppose that it is possible 

to sample A^/,/,rf. samples, called particles, (^oi+i.^o:?+L ) according to 

p(^o:f+z. ® \y\4^L \ Th^D ^ empirical estimate of this distribution is given by 

T7Zr=»^jc<o 9<'> (^o:-.^Qo:,.J (22) 

so that a Monte Carlo approximation of the marginal distribution p (dxi, d9t|yi:i+L) 
follows as 



(23) 
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Using this distribution, an estimate of II (ft) for any ft may be obtained as 



\il*Ll 

1 « (24) 

=^Z:./,(-^'^^/•^)- 



This estimate is unbiased and from the strong law of large numbers, 

^ a.5. 

Il,n (/t) N-^^oo (/O- Under additional assumptions, the estimates satisfy a central limit 
theorem. The advantage of the Monte Carlo method is clear. It is easy to estimate II 
(/ 1) for any /u and the rate of convergence of this estimate does not depend on t or the 
dimension of the state space, but only on the number of particles N and the 
characteristics of the function /t. Unfortunately, it is not possible to sample directly 
from the distribution p (dxon+L, d0o:n^L|yi:t+L) at any t, and alternative strategies need to 
be investigated. 



One solution to estimate p (dxo:t+u<10o:t+Liyi;t+L) and 

II (/t) is the well-known Bayesian importance sampling method as disclosed in A. 
Doucet, S.J. Godsill and C. Andrieu, "On sequential Monte Carlo sampling methods for 
Bayesian filtering". Statistics and Computing, 2000, the contents of which are 
incorporated herein by reference. This method assumes the existence of an arbitrary 
importance distribution n (dxon+u d9o:t+L|yo:i+L) which is easily simulated from, and 
whose support contains that of p (dxo:t+udeo:t+L|yi:t+L)- Using this distribution It {fx) 
may be expressed as 



where the importance weight w (xo:t+L,6o:t+L) is given by 



M^O:/+i »t> p, — j (26) 
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The importance weight can normally only be evaluated up to a constant of 
proportionality, since, following from Bayes* rule 

_ Pjyyj.L Wj.l >9(>..£ )Pidx^,.L > d^o^..L ) (27) 

where the normalizing constant p(yi:t+L) can typically not be expressed in closed-form. 

If N i.i.d. samples {x^^^^^OoL^ ) can be simulated according to a distribution 

7c (dxo:t+L,d9o:t+L|yi:t+L), a Montc Carlo estimate of II (ft) in (25) may be obtained as 



^ L.N ~ — 

(28) 



where the normalized importance weights are given by 



(29) 



This method is equivalent to a point mass approximation of the target distribution of the 
form 



ZN (30) 



The perfect simulation case, when n (dxo:t+L, deo:t+L|yi:t+L) = p (dxo:t+L, deoii+Llyia+t), 
corresponds to Tv^^;^^^ = N"' , i = 1 , . . .,N. In practice, the importance distribution will be 
chosen to be as close as possible to the target distribution in a given sense. For finite N, 
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(fx) is biased, since it involves a ratio of estimates, but asymptotically; according to 
the strong law of large numbers, l[ ^ (ft) 77^ II (ft) Under additional assumptions 
a central limit theorem also holds as disclosed in Doucet, Godsill and Andrieu. 

It is possible to reduce the estimation of p (dxt, d9i|yi:t+L) and II (ft) to one of sampling 
from p (dyo:t+Ltyi:i+L), where we recall that yi A {x, ,a . Indeed, 

p (dai,dYo:t+Lly 1 :t+L) = P (dat|yo:i+u y i ct+O X p (dyo:t+L I y 1 rt+O (31) 

where p (datlyort+t^yiit+L) is a Gaussian distribution whose parameters may be computed 
using Kalman filter type techniques. Thus, given an approximation of p (dyo:t+L|yi:i+L)j 
an approximation of p (dxt,d0t|yi:i+L) may straightforwardly be obtained. Defining the 
marginal importance distribution and associated importance weight as 

Oj+L 

p(y o^,*L\yi..*L) (32) 



^(yO:/.iK/*i.)' 



and assuming that a set of samples y distributed according to n {dy ^^^^ \yi j^L ) is 
available, an alternative Bayesian importance sampling estimate of II (/t) follows as 

KMy'^.J (33) 



provided that p (at|yo:t+L,yi:t+L)/i (xt, 9t) can be evaluated in a closed form expression. In 
(33) the marginal normalized importance weights are given by 
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Intuitively, to reach a given precision, 1 1 will need a reduced number of samples 

over / (/, ) , since it only requires samples from the marginal distribution 

^ o^li l^^i f+L ) It can be proved that the variances of the estimates is subsequently 

reduced as disclosed in A. Doucet, J.F.G. de Freitas and N J. Gordon (eds.). Sequential 
Monte Carlo Methods in Practice, Springer-Verlag, June 2000, the contents of which 
are incorporated herein by reference. In the present case, this is important as at each 
time instant the number of parameters is (when assuming that all mixing filters and AR 
processes have the same length), 

• m^L - mL parameters for the mixing filters, where L can be large. 

• m or 1 parameter(s) for the observation noise. 

• nl -J- n parameters for the autoregressive processes. 

• n parameters for the sources. 

It is not clear which integration will allow for the best variance reduction, but, at least in 
terms of search in the parameters space, the integration of the mixing filters and 
autoregressive filters seems preferable. 

Given these results, the subsequent discussion will focus on Bayesian importance 
sampling methods to obtain approximations of p(dY^^^^\y^.,^^ ) and II (ft) using an 

importance distribution of the form ti (rfy |>;, ) . The methods described up to now 
are batch methods. How a sequential method may be obtained is described below. 

The importance distribution at discrete time t may be factorized as 
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The aim is to obtain at any time t an estimate of the distribution p (<iyo:.+L|yi:t+L) and to 
be able to propagate this estimate in time without modifying subsequently the past 
simulated trajectories y^]^^^ . This means that n (dyo:,+L|yin+L) should admit n (dyo:,. 

i+Llyirt-i+t) as marginal distribution. This is possible if the importance distribution is 
restricted to be of the general form 

Such an importance distribution allows recursive evaluation of the importance weights, 
i.e. w (you+l) = w (yo:t-i+L) w,+l, and in this particular case 



"(Yo:,.ih:/*i) 71(7 



The quantity p (dx,tLix,+i.,p,+t) can be computed up to a normalizing constant using a 
one step ahead Kalman filter for the system given by Eq. (19) and p (yt+tlxo i+L.Pt+L) can 
be computed using a one step ahead Kalman filter of the system given by Eq. (17). 

There is an unlimited number of choices for the importance distribution it (dyo:,+L|yi:t+L), 
the only restriction being that its support includes that of p (dyo:,+L|yi:t+L)- Two 
possibilities are considered next. A possible strategy is to be choose at time t+L the 
importance distribution that minimizes the variance of the importance weights given yo:,. 
1 and y,:,. The importance distribution that satisfies this condition is p (dynLlyo t- 
i+L,yi:t+L), with the associated incremental importance weight given by 

^r^iOt Piy,^L\icy..uL>yv.,.L-l) 

= ^ P(y,*L \f <y.,.L ' yu,*i.i )p{dy (y ). ^^^^ 
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Direct sampling from the optimal importance distribution is difficult, and evaluating the 
importance weight is analytically intractable. The aim is thus in general to mimic the 
optimal distribution by means of tractable approximations, typically local linearization 
P (dyttyo:M,yi:i)- Instead, a mixed suboptimal method is described. We propose to 
sample the particles at time t according to two importance distributions 7i| and with 
proportions a and 1 - a such that the importance weights w (y o'^^^ ) have now the form 
(note that it would be possible to draw Ni and N2 randomly according to a Bernoulli 
distribution with parameter a, but this would increase the estimator variance) 



a 



P(y 0:^.L\yVJ.L) 



O-a) , (.) I 



7ti(Ya,*/.h:,.i). 

p(y o-j^Ayvj^t) 



(39) 



which in practice is estimated as w^.'/^.j 



piy iirVi \yu*L ) / " ^ (y a'/li \yv.,.L ) 

(l-a)- 



s;^., p(y i^lt . ) / 7t Ky o:?.i \y,..L) 



The importance distribution iti (dx,+L!xi:t+L-i,Pi:t+uyi:t+L) will be taken to a normal 
distribution centered around zero with variance a^, and 712 (dx,+L|xi:i+L-i,Pi:t+L,yi:t+L) is 

taken to be p(dx^i\\nt:f^{,\,_ .P/^^i - P,?t -3',..^ ) which is a Gaussian distribution 

obtained from a one step ahead Kalman filter for the state space model described in (14) 
with fw"^'?, and OT*^'?, , as values for a^'^ and h^'^ and initial variances P*,+L|t+L and 

/+L|/ + i t + L\/+L 

P**t+L|t+L- The variances are sampled from their prior distributions, and expression (37) is 
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used to compute w^^^^ . Note that other importance distributions are possible, but this 
approach yields good results and seems to preserve diversity of the samples. 

For importance distributions of the form specified by (36) the variances of the 
importance weights can only increase (stochastically) over time, as disclosed by Doucet, 
Godsill and Andrieu and the references therein. It is thus impossible to avoid a 
degeneracy phenomenon. Practically, after a few iterations of the algorithm, all but one 
of the normalized importance weights are very close to zero, and a large computational 
effort is devoted to updating trajectories whose contribution to the final estimate is 
almost zero. For this reason it is of cmcial importance to include selection and 
diversity. This is discussed in more detail hereinafter. 

The purpose of a selection (or resampling) procedure is to discard particles with low 
normalized importance weights and multiply those with high normalized importance 
weights, so as to avoid the degeneracy of the algorithm. A selection procedure 
associates with each particle, say f^] , a number of children N, gN, such that 

I ^1 = ^ » to obtain N new particles y^J . If Nj = 0 then f^] is discarded, otherwise 
it has Ni children at time t+1. After the selection step the normalized importance 
weights for all the particles are reset to thus discarding all information regarding the 
past importance weights. Thus, the normalized importance weight prior to selection in 
the next time step is proportional to (37). These will be denoted as w]'^ , since they do 
not depend on any past values of the normalized importance weights. If the selection 
procedure is performed at each time step, then the approximating distribution before the 



selection step is given by /j^ (rfy 
selection step follows as {dy ^.^ 



) ^ S.-i ^^'^^ -(/) o r and the one after the 



y\j) = ^~'Y!!-y^^'^^ oA^(yjy Systematic sampling 

as disclosed by Doucet, de Freitus and Gordon is chosen for its good variance 
properties. 
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However selection poses another problem. During the resampling stage, any particular 
particle with a high importance weight will be duplicated many times. In particular, 
when L>0, the trajectories are resampled L times from time t +1 to t + L so that very 
few distinct trajectories remain at time t + L. This is the classical problem of depletion 
of samples. As a result the cloud of particles may eventually collapse to a single 
particle. This degeneracy leads to poor approximation of the distributions of interest. 
Several suboptimal methods have been proposed to overcome this problem and 
introduce diversity amongst the particles. Most of these are based on kernel density 
methods (as disclosed by Doucet, Godsill and Andrieu and by N.J. Gordon, D.J. 
Salmond and A.F.M. Smith, "Novel approach to nonlinear/non-Gaussian Bayesian state 
estimation", lEE Proceedings-F, vol. 140, no. 2, pp. 107-1 13, 1993, (the contents of 
which are incorporated herein by reference) which approximate the probability 
distribution using a kernel density estimate based on the current set of particles, and 
sample a new set of distinct particles from it. However, the choice and configuration of 
a specific kernel are not always straightforward. Moreover, these methods introduce 
additional Monte Cario variation. It is shown hereinafter how MCMC methods may be 
combined with sequential importance sampling to introduce diversity amongst the 
samples without increasing the Monte Carlo variation. 

An efficient way of limiting sample depletion consists of simply adding an MCMC step 
to the simulation-based filter/fixed-lag smoother (see Berzuini and Gilks referred to by 
Doucet, Godsill and Andrieu and by CP. Robert and G. Casella, Monte Cario Statistical 
Methods, Springer Verlag, 1999, the contents of which are incorporated herein by 
reference). This introduces diversity amongst the samples and thus drastically reduces 
the problem of depletion of samples. Assume that, at time t + L, the particles Xo:ml 
marginally distributed according to p (dyoit+Uyirt+L)- If a transition kernel K 
(Yo:i+L|dYo:t+L) With invariant distribution p (dYo:i+iJyi:t+L) is applied to each of the 
particles, then the new particles y^'l^^ are still distributed according to the distribution 
of interest- Any of the standard MCMC methods, such as the Metropolis-Hastings 
(MH) algorithm or Gibbs sampler, may be used. However, contrary to classical MCMC 
methods, the transition kernel does not need to be ergodic. Not only does this method 
introduce no additional Monte Cario variation, but it improves the estimates in the sense 
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that it can only reduce the total variation norm of the current distribution of the particles 
with respect to the target distribution. 

Given at time t + L-1, NeN* particles y^]+^_, distributed approximately according to 

P (dyo:t+L-i|yi:t+L-i), the Monte Carlo fixed-lag smoother proceeds as follows at time t + 
L. 

Sequential Importance Sampling Step 

• For i - 1 , . . . .,N, f /;>, ~ 71 (rfy , y,,,, ) and set f^l^ = (f i^,^., )- 

• For i == 1, . ..,,N, compute the normalized importance weights >v^^^\ using (37) and 
(40). 

Selection Step 

• Multiply / discard particles y^]^i^ with respect to high / low normalised importance 
weights to obtain N particles y ^ ^-S- using systematic sampling, 

MCMC Step 

• For i = 1 , . , ..,N, apply to y '^^.l ^ M^ov transition kemel K (y ^'/^^ jrfy ^^/^^ ) with 
invarient distribution p (dyo:t+L|yi:t+L) to obtain N particles y ^'/^^ . 

There is an unlimited number of choices for the MCMC transition kemeL Here a one- 
at-a-time MH algorithm is adopted that updates at time t + L the values of the Markov 
process from time t to t + L. More specifically y |['\ k = t, . . .,t + L, i = 1 . ..,N, is 

sampled according to p (rfy , |y , y,,^, ), with y <2 A (y ,y ,....,y ;<;> ,....,y ]Tl )• " is 
straightforward to verify that this algorithm admits p (dyo:t+L|yi:t+L) as invariant 
distribution. Sampling fh>m p cfy |y I'i >yi:,+t ) can be done efficiently via a backward - 
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forward algorithm of O (L+1) complexity as disclosed in A. Doucet and C. Andrieu, 
"Iterative algorithms for optimal state estimation of jump Markov linear systems" in 
Proc. Conf IEEE ICASSP, 1999 the contents of which are incorporated herein by 
reference. At time t + L it proceeds as summarised below for the i-th particle. 

For k = t + L , t, compute and store P^'^l^m'^^^^^ and P^^'l^ by running the information 

filter defined in (50)-(51) of the Appendix for the two systems (17) and (19). 

Forward Step 

) 

Fork = t, , t4-L. 

• - Sample a proposal Yr ~ q (dYk|y-k,yo:i+L), using the proposal distribution in (43). 

- Perform one step of the Kalman filter in (48)-(49) for the current value y and 

the proposed value y^, and calculate their posterior probabilities using (41) and 
for the two systems (17) and (19). 

- Compute the MH acceptance probability a(y ^^'^ ,y ) , as defined in (44). 

- If (u - Ufo.i]) < a (y ,y J , set y ['^ = y^, otherwise set y J'^ - y . 

The target posterior distribution for each of the MH steps is given by 

p(^,,Pjx_4,p_,,y,^,^J 

)p(x,,p,|x.,.p_J 

)pixt l^-i , Po^.i )p(Pa , p., ) 

« Piyk 1^0:*-, . Po:*-1 ) X lpiyk*U,.L\K'^k*VJ.L^^k^U*L)PiK\yV.,,X^:k . Po:* M^A 

X P(P« |P*-, )P(P |P A )/'(^* |po..i ) Ip{x,^U.L K ' P l^*X:.*L K:* , P 0:* )^«* 
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There is a similarity between the two expressions in yk and Xk. These two terms can be 
computed in an identical way by first using two KaUnan filters for system (17) and (19) 
to obtain p (yk|xo:k-i,Po:ic.i) = Niy, ; y^^,,^S^, ) and p(x Jpo-^, J = N^x^ ; x^^^^^S; ) and 
secondly two information filters to obtain 



jpiy 



and a similar expression for J pix,^^yj^i^\a^ y^t^ij+i)pia,^ [x,^ , Po:*)^^* - The different 

matrices involved are defined as follows. Let P^^^ = R^nlR^^ , where fi* g R"'*''"'* is 
the diagonal matrix containing the «^ < «^ non-zero singular values of 

, and /fjf e i?"*'"* is the matrix containing the columns of R^ corresponding to the 

non-zero singular values, where P^^^ ~ R^UlR^^ is the singular value decomposition of 

^4* matrix is given by 

To sample from the distribution in (41) using a MH step, the proposal distribution is 
here taken to be: 

9(P* |P*-, . P*., )a/'(P*., |P* )KP*|P*-, ) 

Qi^k F-* . Po:,*i .y^+t )ap(^* \m1^t . 'w* . P * '-Ko:* ) 

and 

* hr -* ,3'o:,*i )a^{Pt |P *-i . P*.. k-* . Po.,*i . yo:,*i )• (43) 
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which requires a one step ahead Kalman filter on the system (14). In both cases 

4>jt (4>i-i y^k4ri ( ^ If the current and proposed new values for the 

State of tlie Markov chain are given by yk and y\, respectively, the MH acceptance 
probability follows as 

a (y * , Y I) = min {1, r(y ^ ,y ; )} , (44) 
with the acceptance ratio given by 

At each iteration the computational complexity of the Monte Carlo fixed-lag smoother 
is O ((L + 1)N), and it is necessary to keep in memory the paths of all the trajectories 
from time t to t + L, i.e. {y )!}^^ : / = 1,...., N} as well as the sufficient statistics 

The computational complexity of this algorithm at each iteration is clearly O (N). At 
first glance, it could appear necessary to keep in memory the paths of all the trajectories 
ou ' ' = ^y "'y so that the storage requirements would increase linearly with time. 
In fact, for both the optimal and prior importance distributions, n (yt|yo:t-i>yi:t) and the 
associated importance weights depend on yo:M only via a set of low-dimensional 
sufficient statistics {m^^f^'^ , P^^;^^'^ : i = 1,..,., M}, and only these values need to be kept in 

memory. Thus, the storage requirements are also O (N) and do not increase over time. 
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Appendix: Kalman filter recursions 
The system 

x,^,=A,^xX,+B,,,v,^, (46) 
y,=C.X,+D,w, (47) 



is considered. 



The sequence <;i:t being here assumed known, the Kalman filter equations are the 
following. 

Set mo|o= Qo and Poio =Po> then for t = 1 ,T compute 



Pt\t ~^/|f-J "Pflt-l^t CtPt\t-} 



(47) 



(48) 



where mttt-i = E{xi|yi:t.i,qi:t}, nh|t= E {xt|yii, q^}, yt\t-\ = yi - ytfi-i, Pt|t-i 
cov {xi|yi:i.i, ^i:t}, Pttt = cov { Xi|y i:t,qi:i}, yi(t.i = E {yi|yi:t-i, <;\:t} and 
Si = cov {yi|yi:t-i,^i:t}- The likelihood p (yi|qi:i) is estimated as 



(49) 



The backward information filter proceeds as follows from time t + L to t. 
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andfork = t + L-l,....,l 

p'-l _ aT pr-\ 
^ ^*+i|* + l'"*+l|jfc+l 

^;il'=n'ill,+c;(AZ);)-c, 



CLAIMS: 
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1 . A method of extracting at least one desired signal from a system cx>mprising at 
least one measured contaminated signal and at least one commimication channel via 
which the at least one contaminated signal is measured, comprising modelling the 
system as a state space model in which the at least one desired signal has first 
characteristics and the at least one communication channel has second characteristics 
which are different from the first characteristics. 

2. A method as claimed in claim 1, in which the first characteristics are fixed and 
known. 

3. A method as claimed in claim 1, in which the first characteristics are time- 
varying. 

4. A method as claimed in claim lor 3, in which the second characteristics are 
fixed and known. 

5. A method as claimed in any one of claims 1 to 3, in which the second 
characteristics are time-varying. 

6. A method as claimed in claim 5 when dependent on claim 3, in which the first 
time- varying characteristics vary on average more quickly than the second time- varying 
characteristics. 

7. A method as claimed in any one of the preceding claims, in which the at least 
one communication charmel comprises a plurality of communication channels and the at 
least one contaminated signal comprises a plurality of contaminated signals. 

8. A method as claimed in any one of the preceding claims, in which the at least 
one desired signal comprises a plurality of desired signals. 
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9. A method as claimed in claim 7 or in claim 8 when dependent on claim 7, in 
which the number of communication channels is greater than or equal to the number of 
desired signals. 

10. A method as claimed in claim 8 or in claim 9 when dependent on claim 8, in 
which the at least one contaminated signal comprises a linear combination of time- 
delayed versions of at least some of the desired signals. 

11. A method as claimed in any one of the preceding claims, in which the at. least 
one contaminated signal comprises the at least one desired signal contaminated with 
noise. 

12. A method as claimed in any one of the preceding claims, in which the at least 
one channel comprises a plurality of signal propagation paths of different lengths. 

13. A method as claimed in any one of the preceding claims, in which the at least 
one desired signal comprises a sound signal, 

14. A method as claimed in claim 13, in which the least one sound signal comprises 
speech. 

15. A method as claimed in claim 13 or 14 when dependent on claim 8, in which the 
contaminated signals are measured by spatially sampling a sound field. 

16. A method as claimed in any one of the preceding claims, in which the at least 
one desired signal is modelled as a time-varying autoregression. 

1 7. A method as claimed in any one of claims I to 1 5, in which the at least one 
desired signal is modelled as a moving average model. 

18. A method as claimed in any one of claims 1 to 15, in which the at lezist one 
desired signal is modelled as a non-linear time- varying model. 



37 

19, A method as claimed in any one of the preceding claims, in which the at least 
one communication channel is modelled as a finite impulse response model, 

20. A method as claimed in any one of claims 1 to 1 8, in which the at least one 
communication channel is modelled as an infinite impulse response model. 

2L A method as claimed in any one of claims 1 to 18, in which the at least one 
communication channel is modelled as a non-linear time-varying model. 

22. A method as claimed in any one of the preceding claims, in which the state 
space model has at least one parameter which is modelled using a probability model. 

23. A method as claimed in claim 22, in which the at least one desired signal is 
extracted by a Bayesian inversion, 

24. A method as claimed in claim 23, in which the Bayesian inversion is performed 
by a sequential Monte Carlo method or particle filter. 

25. A program for controlling a computer to perform a method as claimed in any 
one of the preceding claims. 

26. A carrier containing a program as claimed in claim 25. 



27. A computer programmed by a program as claimed in claim 25. 
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